MODULE SUSDP_DERIV_CTL_MOD
CONTAINS
SUBROUTINE SUSDP_DERIV_CTL(KIDIA, KFDIA, PSLT, KLEVS3D, PSSDP2, PSSDP3, YDVEG, YDSOIL, YDAGS, YDCST)

USE PARKIND1, ONLY : JPIM, JPRB, JPRD
USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_VEG , ONLY : TVEG
USE YOS_SOIL, ONLY : TSOIL
USE YOS_AGS , ONLY : TAGS
USE YOS_CST , ONLY : TCST
USE YOMSURF_SSDP_MOD
USE COTWO_MOD

! (C) Copyright 2025- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!     ---------------------------------------------------------------------------
!**   *SUSDP_DERIV* - COMPUTE DERIVED SDP PARAMETERS 

!      PURPOSE
!      -------
!      COMPUTE THE DERIVED SURFACE SPATIALLY DISTRIBUTED PARAMETERS THAT CAN BE
!      CALIBRATED THROUGH THE READING OF AN EXTERNAL FILE 
!
!      INTERFACE
!      ---------
!      *SUSDP_DERIV* IS CALLED BY *SUINIF1S* AND *SUGRIDG*

!      INPUT PARAMETERS (REAL):
!      *PSLT*        SOIL TYPE

!      INPUT PARAMETERS (INTEGER):
!      *KIDIA*       START POINT
!      *KFDIA*       END POINT
!      *KLEVS3D*     NUMBER OF SOIL VERTICAL LEVELS

!      OUTPUT PARAMETERS (REAL):
!      *PSSDP2*      2D OBJECT WITH SURFACE SPATIALLY DISTRIBUTED PARAMETERS
!      *PSSDP3*      3D OBJECT WITH SURFACE SPATIALLY DISTRIBUTED PARAMETERS
!     

!      REFERENCE
!      ---------

!      Original:   Iria Ayan-Miguez (BSC) September 2023

!     --------------------------------------------------------------------------


IMPLICIT NONE

INTEGER(KIND=JPIM), INTENT(IN)  :: KIDIA
INTEGER(KIND=JPIM), INTENT(IN)  :: KFDIA
REAL(KIND=JPRB),    INTENT(IN)  :: PSLT(:)
REAL(KIND=JPRB),    INTENT(OUT) :: PSSDP2(:,:)
REAL(KIND=JPRB),    INTENT(OUT) :: PSSDP3(:,:,:)
INTEGER(KIND=JPIM), INTENT(IN)  :: KLEVS3D
TYPE(TVEG),         INTENT(IN)  :: YDVEG
TYPE(TSOIL),        INTENT(IN)  :: YDSOIL
TYPE(TAGS),         INTENT(IN)  :: YDAGS
TYPE(TCST),         INTENT(IN)  :: YDCST

INTEGER(KIND=JPIM) :: JK, JL
REAL(KIND=JPRB) :: ZCONV, ZPSICAP, ZPSIPWP, ZFAC, ZRHOD
REAL(KIND=JPRB) :: ZQ, ZLAMBDASM, ZEPSILON, ZWMAX, ZWMIN, ZWFAC_S
REAL(KIND=JPRB) :: ZAG(1),ZRD(1),ZGS(1), ZF2, ZFZEROSTARH(1), ZRVFZEROST(1)
REAL(KIND=JPRB) :: ZAMMAXTL(1), ZAMMAXTH(1), ZGAMMTL(1),ZGAMMTH(1), ZGMESTL(1), ZGMESTH(1)
REAL(KIND=JPRB) :: ZRCO2(1),ZRDSPOPT(1),ZRIAOPT(1),ZRVANMAX(1),ZRVGC(1),ZRVDMAX(1),ZRVEPSO(1)
REAL(KIND=JPRD) :: ZSEMAX, ZDMAX_D, ZSEMIN, ZDMIN_D

! local variables for COTWO call for one grid
INTEGER(KIND=JPIM) :: ZIDIA=1_JPIM, ZFDIA=1_JPIM, ZLON=1_JPIM
LOGICAL :: ZLDLAND(1)=.TRUE.


REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

IF (LHOOK) CALL DR_HOOK('SUSDP_DERIV_CTL_MOD:SUSDP_DERIV_CTL',0,ZHOOK_HANDLE)
ASSOCIATE(RVHSTRH2D=>PSSDP2(:,SSDP2D_ID%NRVHSTRH2D), RVHSTRL2D=>PSSDP2(:,SSDP2D_ID%NRVHSTRL2D),& 
    & RVSEFOLDH2D=>PSSDP2(:,SSDP2D_ID%NRVSEFOLDH2D), RVSEFOLDL2D=>PSSDP2(:,SSDP2D_ID%NRVSEFOLDL2D),&   
    & RMVGALPHA3D=>PSSDP3(:,:,SSDP3D_ID%NRMVGALPHA3D), RNFACM3D=>PSSDP3(:,:,SSDP3D_ID%NRNFACM3D),&
    & RWRESTM3D=>PSSDP3(:,:,SSDP3D_ID%NRWRESTM3D), RWSATM3D=>PSSDP3(:,:,SSDP3D_ID%NRWSATM3D),&
    & RVAMMAXH2D=>PSSDP2(:,SSDP2D_ID%NRVAMMAXH2D), RVAMMAXL2D=>PSSDP2(:,SSDP2D_ID%NRVAMMAXL2D),&
    & RVCEH2D=>PSSDP2(:,SSDP2D_ID%NRVCEH2D), RVCEL2D=>PSSDP2(:,SSDP2D_ID%NRVCEL2D),&
    & RVCFH2D=>PSSDP2(:,SSDP2D_ID%NRVCFH2D), RVCFL2D=>PSSDP2(:,SSDP2D_ID%NRVCFL2D),&
    & RVCNAH2D=>PSSDP2(:,SSDP2D_ID%NRVCNAH2D), RVCNAL2D=>PSSDP2(:,SSDP2D_ID%NRVCNAL2D),&
    & RVDMAXH2D=>PSSDP2(:,SSDP2D_ID%NRVDMAXH2D), RVDMAXL2D=>PSSDP2(:,SSDP2D_ID%NRVDMAXL2D),&
    & RVEPSOH2D=>PSSDP2(:,SSDP2D_ID%NRVEPSOH2D), RVEPSOL2D=>PSSDP2(:,SSDP2D_ID%NRVEPSOL2D),&
    & RVFZEROSTH2D=>PSSDP2(:,SSDP2D_ID%NRVFZEROSTH2D), RVFZEROSTL2D=>PSSDP2(:,SSDP2D_ID%NRVFZEROSTL2D),&
    & RVGAMMH2D=>PSSDP2(:,SSDP2D_ID%NRVGAMMH2D), RVGAMML2D=>PSSDP2(:,SSDP2D_ID%NRVGAMML2D),&
    & RVGCH2D=>PSSDP2(:,SSDP2D_ID%NRVGCH2D), RVGCL2D=>PSSDP2(:,SSDP2D_ID%NRVGCL2D),&
    & RVGMESH2D=>PSSDP2(:,SSDP2D_ID%NRVGMESH2D), RVGMESL2D=>PSSDP2(:,SSDP2D_ID%NRVGMESL2D),&
    & RVQDAMMAXH2D=>PSSDP2(:,SSDP2D_ID%NRVQDAMMAXH2D), RVQDAMMAXL2D=>PSSDP2(:,SSDP2D_ID%NRVQDAMMAXL2D),&   
    & RVQDGAMMH2D=>PSSDP2(:,SSDP2D_ID%NRVQDGAMMH2D), RVQDGAMML2D=>PSSDP2(:,SSDP2D_ID%NRVQDGAMML2D),&
    & RVQDGMESH2D=>PSSDP2(:,SSDP2D_ID%NRVQDGMESH2D), RVQDGMESL2D=>PSSDP2(:,SSDP2D_ID%NRVQDGMESL2D),&
    & RVT1AMMAXH2D=>PSSDP2(:,SSDP2D_ID%NRVT1AMMAXH2D), RVT1AMMAXL2D=>PSSDP2(:,SSDP2D_ID%NRVT1AMMAXL2D),& 
    & RVT1GMESH2D=>PSSDP2(:,SSDP2D_ID%NRVT1GMESH2D), RVT1GMESL2D=>PSSDP2(:,SSDP2D_ID%NRVT1GMESL2D),&
    & RVT2AMMAXH2D=>PSSDP2(:,SSDP2D_ID%NRVT2AMMAXH2D), RVT2AMMAXL2D=>PSSDP2(:,SSDP2D_ID%NRVT2AMMAXL2D),& 
    & RVT2GMESH2D=>PSSDP2(:,SSDP2D_ID%NRVT2GMESH2D), RVT2GMESL2D=>PSSDP2(:,SSDP2D_ID%NRVT2GMESL2D),&
    & RVTOPTH2D=>PSSDP2(:,SSDP2D_ID%NRVTOPTH2D), RVTOPTL2D=>PSSDP2(:,SSDP2D_ID%NRVTOPTL2D),&
    & RVANMAXH2D=>PSSDP2(:,SSDP2D_ID%NRVANMAXH2D), RVANMAXL2D=>PSSDP2(:,SSDP2D_ID%NRVANMAXL2D),&
    & RVBSLAI_NITROH2D=>PSSDP2(:,SSDP2D_ID%NRVBSLAI_NITROH2D), RVBSLAI_NITROL2D=>PSSDP2(:,SSDP2D_ID%NRVBSLAI_NITROL2D),&
    & RLAMBDADRYM3D=>PSSDP3(:,:,SSDP3D_ID%NRLAMBDADRYM3D), RLAMSAT1M3D=>PSSDP3(:,:,SSDP3D_ID%NRLAMSAT1M3D), &
    & RRCSOILM3D=>PSSDP3(:,:,SSDP3D_ID%NRRCSOILM3D), RWCAPM3D=>PSSDP3(:,:,SSDP3D_ID%NRWCAPM3D),&
    & RWPWPM3D=>PSSDP3(:,:,SSDP3D_ID%NRWPWPM3D), RCGDRYM3D=>PSSDP3(:,:,SSDP3D_ID%NRCGDRYM3D), &
    & RMFACM3D=>PSSDP3(:,:,SSDP3D_ID%NRMFACM3D), RLAMBDAM3D=>PSSDP3(:,:,SSDP3D_ID%NRLAMBDAM3D), &
    & RDMAXM3D=>PSSDP3(:,:,SSDP3D_ID%NRDMAXM3D), RDMINM3D=>PSSDP3(:,:,SSDP3D_ID%NRDMINM3D), &
    & RWCONSM3D=>PSSDP3(:,:,SSDP3D_ID%NRWCONSM3D), &
    & RETV=>YDCST%RETV, LEVGEN=>YDSOIL%LEVGEN,RVGBARCAP=>YDSOIL%RVGBARCAP,RBARCAP=>YDSOIL%RBARCAP,&
    & RGH2O=>YDSOIL%RGH2O, RCO2=>YDAGS%RCO2, RDSPOPT=>YDAGS%RDSPOPT, RIAOPT=>YDAGS%RIAOPT,&
    & RAW=>YDAGS%RAW, RBW=>YDAGS%RBW, RDAY=>YDCST%RDAY, RRHOSM=>YDSOIL%RRHOSM, &
    & RLAMBDAQ=>YDSOIL%RLAMBDAQ, RLAMBDAO=>YDSOIL%RLAMBDAO, RBARPWP=>YDSOIL%RBARPWP,&
    & RLAMBDRYMA=>YDSOIL%RLAMBDRYMA, RLAMBDRYMB=>YDSOIL%RLAMBDRYMB, RLAMBDRYMC=>YDSOIL%RLAMBDRYMC)
 

! VEGETATION PARAMETERS
ZCONV=1013.25_JPRB*(RETV+1)

RVHSTRH2D(:)=ZCONV*RVHSTRH2D(:)
RVHSTRL2D(:)=ZCONV*RVHSTRL2D(:)

! AGS PARAMETERS

IF (YDVEG%LECTESSEL) THEN
  
  ! RVSEFOLD e-folding time for senescence (s)
  RVSEFOLDL2D(:)=RDAY*RVSEFOLDL2D(:)
  RVSEFOLDH2D(:)=RDAY*RVSEFOLDH2D(:)

  ! from nitrogen decline theory
  RVBSLAI_NITROH2D(:) = 1._JPRB/(RVCEH2D(:)*RVCNAH2D(:)+RVCFH2D(:))
  RVBSLAI_NITROL2D(:) = 1._JPRB/(RVCEL2D(:)*RVCNAL2D(:)+RVCFL2D(:))

   
  ! Calculate RVANMAX (Compute optimum (maximum) net CO2 assimilation (An,max))
  ! Initialize local variables
  ZAG(1)=0.0_JPRB   ! Agross (diagn)
  ZRD(1)=0.0_JPRB   ! Rd (diagn)
  ZGS(1)=0.0_JPRB   ! leaf conductance to H20 (m s-1)
  ZAMMAXTL(1)=0.0_JPRB ! maximum photosynthetic capacity low veg (kgCO2 kgAir-1 m s-1)
  ZAMMAXTH(1)=0.0_JPRB ! maximum photosynthetic capacity high veg (kgCO2 kgAir-1 m s-1)
  ZGAMMTL(1)=0.0_JPRB ! CO2 compensation point at T=ZTOPT low veg (kgCO2 kgAir-1)
  ZGAMMTH(1)=0.0_JPRB ! CO2 compensation point at T=ZTOPT high veg (kgCO2 kgAir-1)
  ZGMESTL(1)=0.0_JPRB ! mesophyll conductance at T=ZTOPT low veg (m s-1)
  ZGMESTH(1)=0.0_JPRB ! mesophyll conductance at T=ZTOPT high veg (m s-1)
  ZFZEROSTARH(1)=0.0_JPRB ! f0 unstressed high veg
  ZRCO2(1)    = RCO2
  ZRDSPOPT(1) = RDSPOPT
  ZRIAOPT(1)  = RIAOPT

  DO JL=1,KFDIA
  ! compute temperature responses:
  ! low vegetation
    ZGAMMTL(1)  = RVGAMML2D(JL)*(RVQDGAMML2D(JL)**(0.1_JPRB*(RVTOPTL2D(JL)-25.0_JPRB)))
    ZAMMAXTL(1) = ( RVAMMAXL2D(JL)*RVQDAMMAXL2D(JL)**(0.1_JPRB*(RVTOPTL2D(JL)-25.0_JPRB))) &
               &         /((1.0_JPRB+EXP(0.3_JPRB*(RVT1AMMAXL2D(JL)-RVTOPTL2D(JL))))*          &
               &           (1.0_JPRB+EXP(0.3_JPRB*(RVTOPTL2D(JL)-RVT2AMMAXL2D(JL)))))
    ZGMESTL(1)  = ( RVGMESL2D(JL)*RVQDGMESL2D(JL)**(0.1_JPRB*(RVTOPTL2D(JL)-25.0_JPRB)))  &
               &         /((1.0_JPRB+EXP(0.3_JPRB*(RVT1GMESL2D(JL)-RVTOPTL2D(JL))))*           &
               &           (1.0_JPRB+EXP(0.3_JPRB*(RVTOPTL2D(JL)-RVT2GMESL2D(JL)))))
  ! high vegetation
    ZGAMMTH(1)  = RVGAMMH2D(JL)*(RVQDGAMMH2D(JL)**(0.1_JPRB*(RVTOPTH2D(JL)-25.0_JPRB)))
    ZAMMAXTH(1) = ( RVAMMAXH2D(JL)*RVQDAMMAXH2D(JL)**(0.1_JPRB*(RVTOPTH2D(JL)-25.0_JPRB))) &
               &         /((1.0_JPRB+EXP(0.3_JPRB*(RVT1AMMAXH2D(JL)-RVTOPTH2D(JL))))*          &
               &           (1.0_JPRB+EXP(0.3_JPRB*(RVTOPTH2D(JL)-RVT2AMMAXH2D(JL)))))
    ZGMESTH(1)  = ( RVGMESH2D(JL)*RVQDGMESH2D(JL)**(0.1_JPRB*(RVTOPTH2D(JL)-25.0_JPRB)))  &
               &         /((1.0_JPRB+EXP(0.3_JPRB*(RVT1GMESH2D(JL)-RVTOPTH2D(JL))))*           &
               &           (1.0_JPRB+EXP(0.3_JPRB*(RVTOPTH2D(JL)-RVT2GMESH2D(JL)))))

  ! Add soil moisture stress effect to leaf conductance:
    ZF2=1.0_JPRB ! optimum soil moisture stress factor (-) 
    ZGMESTL(1)=ZGMESTL(1)*ZF2
    ZGMESTH(1)=ZGMESTH(1)*ZF2
    ZFZEROSTARH(1)=(RAW-LOG(RVGMESH2D(JL)*1000._JPRB))/RBW

  ! compute maximum/initial/optimum net assimilation of CO2:  
    ZRVANMAX(1)   = 0.0_JPRB
    ZRVGC(1)      = RVGCL2D(JL)
    ZRVDMAX(1)    = RVDMAXL2D(JL)
    ZRVFZEROST(1) = RVFZEROSTL2D(JL)
    ZRVEPSO(1)    = RVEPSOL2D(JL)
    CALL COTWO(ZIDIA,ZFDIA,ZLON,ZLDLAND, &
            & ZRVANMAX, &
            & ZAG,ZRD,ZGS, &
            & ZRVGC,ZRCO2, &
            & ZRDSPOPT,ZRVDMAX,ZRIAOPT, &
            & ZGAMMTL,YDAGS, &
            & ZRVFZEROST,ZGMESTL,ZRVEPSO,ZAMMAXTL)
    RVANMAXL2D(JL) = ZRVANMAX(1)

    ZRVANMAX(1)   = 0.0_JPRB
    ZRVGC(1)      = RVGCH2D(JL)
    ZRVDMAX(1)    = RVDMAXH2D(JL)
    ZRVEPSO(1)    = RVEPSOH2D(JL)
    CALL COTWO(ZIDIA,ZFDIA,ZLON,ZLDLAND, &
            & ZRVANMAX, &
            & ZAG,ZRD,ZGS, &
            & ZRVGC,ZRCO2, &
            & ZRDSPOPT,ZRVDMAX,ZRIAOPT, &
            & ZGAMMTH,YDAGS, &
            & ZFZEROSTARH,ZGMESTH,ZRVEPSO,ZAMMAXTH)
    RVANMAXH2D(JL) = ZRVANMAX(1)

  ENDDO
ENDIF

! SOIL PARAMETERS

ZFAC=100._JPRB/9.8_JPRB  ! Bit optimized version of the original
ZPSIPWP=RBARPWP*ZFAC   ! Classic permanent wilting point

IF (LEVGEN) THEN
  ZPSICAP=RVGBARCAP*ZFAC  ! Value valid for sandy soil in the literature---> larger AWC 
                           ! favourably compared to FC obs for medium soils (Ukraine/Russia)
                           ! produce invariant equilibrium after rescaling from TESSEL
                           ! enlarge the range of action for SM analysis (active if PWP<W<CAP)
ELSE
  ZPSICAP=RBARCAP*ZFAC  ! Value mostly used in literature (Hillel, 1998)
ENDIF


! Constants for Thermal Diffusivity
! Values taken from Peters-Liddard et al. 1998
! common to CH and MV

ZQ=0.4_JPRB
ZLAMBDASM=(RLAMBDAQ**ZQ)*(RLAMBDAO**(1.0_JPRB-ZQ))
ZEPSILON=100._JPRB*EPSILON(ZEPSILON)

DO JL=KIDIA, KFDIA
  IF (PSLT(JL) >= 1) THEN    

    RWCAPM3D(JL,:) = RWRESTM3D(JL,:)+(RWSATM3D(JL,:)-RWRESTM3D(JL,:)) &
      &       *(1._JPRB/(1._JPRB+((ABS(RMVGALPHA3D(JL,:)*ZPSICAP)) &
      &       **RNFACM3D(JL,:))))**(1._JPRB-(1._JPRB/RNFACM3D(JL,:)))

    RWPWPM3D(JL,:) = RWRESTM3D(JL,:)+(RWSATM3D(JL,:)-RWRESTM3D(JL,:)) &
      &       *(1._JPRB/(1._JPRB+((ABS(RMVGALPHA3D(JL,:)*ZPSIPWP)) &
      &       **RNFACM3D(JL,:))))**(1._JPRB-(1._JPRB/RNFACM3D(JL,:)))
   
    RMFACM3D(JL,:) = 1._JPRB-(1._JPRB/RNFACM3D(JL,:))

    DO JK=1,KLEVS3D
      ZWFAC_S=SIGN(MAX(ABS(RWSATM3D(JL,JK)-RWRESTM3D(JL,JK)),ZEPSILON),(RWSATM3D(JL,JK)-RWRESTM3D(JL,JK)))
      ZWMAX=.999_JPRB*RWSATM3D(JL,JK)
      ZSEMAX =(ZWMAX-RWRESTM3D(JL,JK))/ZWFAC_S
      ZDMAX_D=(ZSEMAX**(RLAMBDAM3D(JL,JK)-(1._JPRB/RMFACM3D(JL,JK)))) &
        &    *(((1._JPRB-(ZSEMAX**(1._JPRB/RMFACM3D(JL,JK))))**(RMFACM3D(JL,JK))) &
        &    + ((1._JPRB-(ZSEMAX**(1._JPRB/RMFACM3D(JL,JK))))**(-RMFACM3D(JL,JK)))-2._JPRB)

      ZWMIN=1.001_JPRB*RWRESTM3D(JL,JK)   
      ZSEMIN=(ZWMIN-RWRESTM3D(JL,JK))/ZWFAC_S 
      ZDMIN_D= (ZSEMIN**(RLAMBDAM3D(JL,JK)-(1._JPRB/RMFACM3D(JL,JK)))) &
        &         *(((1._JPRB-(ZSEMIN**(1._JPRB/RMFACM3D(JL,JK))))**(RMFACM3D(JL,JK))) &
        &         + ((1._JPRB-(ZSEMIN**(1._JPRB/RMFACM3D(JL,JK))))**(-RMFACM3D(JL,JK)))-2._JPRB)
 
      RDMAXM3D(JL,JK) = (((1._JPRB-RMFACM3D(JL,JK))*RWCONSM3D(JL,JK))/ &
        &      (RMVGALPHA3D(JL,JK)*RMFACM3D(JL,JK)*ZWFAC_S)) * REAL(ZDMAX_D,JPRB)
      RDMINM3D(JL,JK) = (((1._JPRB-RMFACM3D(JL,JK))*RWCONSM3D(JL,JK))/ &
        &      (RMVGALPHA3D(JL,JK)*RMFACM3D(JL,JK)*ZWFAC_S)) * REAL(ZDMIN_D,JPRB)
    ENDDO


  ELSE
    RWCAPM3D(JL,:)   = 0.0_JPRB
    RWPWPM3D(JL,:)   = 0.0_JPRB
    RMFACM3D(JL,:)   = 0.0_JPRB
    RDMAXM3D(JL,:)   = 0.0_JPRB
    RDMINM3D(JL,:)   = 0.0_JPRB
  ENDIF


! CLAPP-HORNBERGER

  DO JK=1,KLEVS3D
    ZRHOD = (1.0_JPRB-RWSATM3D(JL,JK))*RRHOSM
    RLAMBDADRYM3D(JL,JK) = (RLAMBDRYMA*ZRHOD+RLAMBDRYMB)/(RRHOSM-RLAMBDRYMC*ZRHOD)  
  ENDDO
  
ENDDO
RLAMSAT1M3D(:,:) = ZLAMBDASM**(1.0_JPRB-RWSATM3D(:,:))

! Soil Heat Capacity
RRCSOILM3D(:,:) = (1.0_JPRB-RWSATM3D(:,:))*RCGDRYM3D(:,:)+RWCAPM3D(:,:)*RGH2O     

END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('SUSDP_DERIV_CTL_MOD:SUSDP_DERIV_CTL',1,ZHOOK_HANDLE)
END SUBROUTINE SUSDP_DERIV_CTL
END MODULE SUSDP_DERIV_CTL_MOD
